{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8b892c8",
   "metadata": {},
   "source": [
    "Printed and electronic copies of *Modeling and Simulation in Python* are available from [No Starch Press](https://nostarch.com/modeling-and-simulation-python) and [Bookshop.org](https://bookshop.org/p/books/modeling-and-simulation-in-python-allen-b-downey/17836697?ean=9781718502161) and [Amazon](https://amzn.to/3y9UxNb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "combined-semiconductor",
   "metadata": {},
   "source": [
    "# Modeling Vaccination"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-table",
   "metadata": {
    "tags": []
   },
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "formal-context",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/' +\n",
    "         'modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "progressive-typing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "breathing-hamilton",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/chap11.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "growing-sperm",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import code from previous notebooks\n",
    "\n",
    "from chap11 import make_system\n",
    "from chap11 import update_func\n",
    "from chap11 import run_simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "identical-steam",
   "metadata": {},
   "source": [
    "In the previous chapter I presented the Kermack-McKendrick (KM) model of infectious disease and used it to model the Freshman Plague at Olin. In this chapter we'll consider metrics intended to quantify the effects of the disease and interventions intended to reduce those effects.\n",
    "\n",
    "We'll use some of the functions from the previous chapter: `make_system`, `update_func`, and the last version of `run_simulation`, which returns the results in a `DataFrame` object."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "complex-renewal",
   "metadata": {},
   "source": [
    "## Immunization\n",
    "\n",
    "Models like this are useful for testing \"what if?\" scenarios. As an\n",
    "example, we'll consider the effect of immunization.\n",
    "\n",
    "Suppose there is a vaccine that causes a student to become immune to the Freshman Plague without being infected. How might you modify the model to capture this effect?\n",
    "\n",
    "One option is to treat immunization as a shortcut from susceptible to\n",
    "recovered without going through infectious. We can implement this\n",
    "feature like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "recent-cooper",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_immunization(system, fraction):\n",
    "    system.init.s -= fraction\n",
    "    system.init.r += fraction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "arranged-screening",
   "metadata": {},
   "source": [
    "`add_immunization` moves the given fraction of the population from `S`\n",
    "to `R`.\n",
    "\n",
    "As a basis for comparison, I'll run the model with the same parameters as in the previous chapter, with no immunization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "found-learning",
   "metadata": {},
   "outputs": [],
   "source": [
    "tc = 3             # time between contacts in days \n",
    "tr = 4             # recovery time in days\n",
    "\n",
    "beta = 1 / tc      # contact rate in per day\n",
    "gamma = 1 / tr     # recovery rate in per day\n",
    "\n",
    "system = make_system(beta, gamma)\n",
    "results = run_simulation(system, update_func)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unsigned-joseph",
   "metadata": {},
   "source": [
    "Now let's see what happens if 10% of students are immune.\n",
    "I'll make another `System` object with the same parameters, then run `add_immunization` to modify the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "enormous-abortion",
   "metadata": {},
   "outputs": [],
   "source": [
    "system2 = make_system(beta, gamma)\n",
    "add_immunization(system2, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "subject-ideal",
   "metadata": {},
   "source": [
    "Now we can run the simulation like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "funny-copper",
   "metadata": {},
   "outputs": [],
   "source": [
    "results2 = run_simulation(system2, update_func)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bronze-techno",
   "metadata": {},
   "source": [
    "The following figure shows `s` as a function of time, with and\n",
    "without immunization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "divided-biotechnology",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.s.plot(style='--', label='No immunization')\n",
    "results2.s.plot(label='10% immunization')\n",
    "\n",
    "decorate(xlabel='Time (days)',\n",
    "         ylabel='Fraction of population')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "passive-dance",
   "metadata": {},
   "source": [
    "With immunization, there is a smaller change in `s`; that is, fewer people are infected.\n",
    "In the next section we'll compute this change and use it to quantify the effect of immunization."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "postal-cemetery",
   "metadata": {},
   "source": [
    "## Metrics\n",
    "\n",
    "When we plot a time series, we get a view of everything that happened\n",
    "when the model ran, but often we want to boil it down to a few numbers\n",
    "that summarize the outcome. These summary statistics are called\n",
    "*metrics*.\n",
    "\n",
    "In the KM model, we might want to know the time until the peak of the\n",
    "outbreak, the number of people who are sick at the peak, the number of\n",
    "students who will still be sick at the end of the semester, or the total number of students who get sick at any point.\n",
    "\n",
    "As an example, I will focus on the last one --- the total number of sick students --- and we will consider interventions intended to minimize it.\n",
    "\n",
    "We can get the total number of infections by computing the difference in `s` at the beginning and the end of the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "synthetic-element",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def calc_total_infected(results, system):\n",
    "    s_0 = results.s[0]\n",
    "    s_end = results.s[system.t_end]\n",
    "    return s_0 - s_end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "parallel-pipeline",
   "metadata": {},
   "source": [
    "And here are the results from the two simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "recovered-picnic",
   "metadata": {},
   "outputs": [],
   "source": [
    "calc_total_infected(results, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "american-transfer",
   "metadata": {},
   "outputs": [],
   "source": [
    "calc_total_infected(results2, system2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adverse-trance",
   "metadata": {},
   "source": [
    "Without immunization, almost 47% of the population gets infected at some point. With 10% immunization, only 31% get infected. That's pretty good."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eight-maximum",
   "metadata": {},
   "source": [
    "## Sweeping Immunization\n",
    "\n",
    "Now let's see what happens if we administer more vaccines. This\n",
    "following function sweeps a range of immunization rates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "progressive-architect",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sweep_immunity(fraction_array):\n",
    "    sweep = SweepSeries()\n",
    "\n",
    "    for fraction in fraction_array:\n",
    "        system = make_system(beta, gamma)\n",
    "        add_immunization(system, fraction)\n",
    "        results = run_simulation(system, update_func)\n",
    "        sweep[fraction] = calc_total_infected(results, system)\n",
    "\n",
    "    return sweep"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "timely-industry",
   "metadata": {},
   "source": [
    "The parameter of `sweep_immunity` is an array of immunization rates. The result is a `SweepSeries` object that maps from each immunization rate to the resulting fraction of students ever infected.\n",
    "\n",
    "We can call it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "measured-pavilion",
   "metadata": {},
   "outputs": [],
   "source": [
    "fraction_array = linspace(0, 1, 21)\n",
    "infected_sweep = sweep_immunity(fraction_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "indie-seeker",
   "metadata": {},
   "source": [
    "The following figure plots the `SweepSeries`. Notice that the $x$-axis is the immunization rate, not time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "interior-humanitarian",
   "metadata": {},
   "outputs": [],
   "source": [
    "infected_sweep.plot(color='C2')\n",
    "\n",
    "decorate(xlabel='Fraction immunized',\n",
    "         ylabel='Total fraction infected',\n",
    "         title='Fraction infected vs. immunization rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "turkish-mumbai",
   "metadata": {},
   "source": [
    "As the immunization rate increases, the number of infections drops\n",
    "steeply. If 40% of the students are immunized, fewer than 4% get sick.\n",
    "That's because immunization has two effects: it protects the people who get immunized (of course) but it also protects the rest of the\n",
    "population.\n",
    "\n",
    "Reducing the number of \"susceptibles\" and increasing the number of\n",
    "\"resistants\" makes it harder for the disease to spread, because some\n",
    "fraction of contacts are wasted on people who cannot be infected. This\n",
    "phenomenon is called *herd immunity*, and it is an important element\n",
    "of public health (see <http://modsimpy.com/herd>)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "french-spouse",
   "metadata": {},
   "source": [
    "The steepness of the curve is a blessing and a curse. It's a blessing\n",
    "because it means we don't have to immunize everyone, and vaccines can\n",
    "protect the \"herd\" even if they are not 100% effective.\n",
    "\n",
    "But it's a curse because a small decrease in immunization can cause a\n",
    "big increase in infections. In this example, if we drop from 80%\n",
    "immunization to 60%, that might not be too bad. But if we drop from 40% to 20%, that would trigger a major outbreak, affecting more than 15% of the population. For a serious disease like measles, just to name one, that would be a public health catastrophe."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "amino-excess",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In general, models are used to predict, explain, and design.\n",
    "In this chapter, we use an SIR model to predict the effect of immunization and to explain the phenomenon of herd immunity.\n",
    "\n",
    "In the repository for this book, you will find a file called *plague.ipynb* that uses this model for design, that is, for making public health decisions intended to achieve a goal.\n",
    "\n",
    "In the next chapter, we'll explore the SIR model further by sweeping the parameters.\n",
    "\n",
    "But first you might want to work on this exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "institutional-memory",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "You can access the notebooks at <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "drawn-hindu",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    " Suppose we have the option to quarantine infected students.  For example, a student who feels ill might be moved to an infirmary or a private dorm room until they are no longer infectious.\n",
    "\n",
    "How might you incorporate the effect of quarantine in the SIR model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "impressive-librarian",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "assumed-license",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "intended-premium",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "limiting-interest",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "undefined-treasury",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
